%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "Cyclical skill-biased technological change"
% Almut Balleer and Thijs van Rens
% Estimation of shocks from simulated data (Lindquist model)
% This version: May 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% data

% choose rho in the production function
for prodfnc = 1:2
    prodfnc;

    switch prodfnc
        case 1
            load('VARdata_067');
        case 2
            load('VARdata_117')
        case 3
            load('VARdata_167')
        case 4
            load('VARdata_217')
        case 5
            load('VARdata_267')
        case 6
            load('VARdata_317')
        case 7
            load('VARdata_500')
    end

    % simulated data size is 110 x 9 x 1000
    % variables: average hours, output per hour, investment price, skill
    % premium, wage unskilled, wage skilled, relative supply (hours), total
    % hours unskilled, total hours skilled workers
    % in levels

    HORIZON_ident = 80;    % identification horizon
    simsize = size(VARdata,3);
    num_vars = 4;
    do_trend = 0;
    Respmedall = zeros(num_vars^2,HORIZON_ident,simsize);
    Resplowall = zeros(num_vars^2,HORIZON_ident,simsize);
    Resphighall = zeros(num_vars^2,HORIZON_ident,simsize);
    tic
    for it=1:simsize
        display(it);
        % assign data
        % no relative supply
        wage_premium = log(VARdata(:,4,it));
        price_series = log(VARdata(:,3,it));
        prod_series = log(VARdata(:,2,it));
        hours_series = log(VARdata(:,1,it));

        d_price = diff(price_series)*100;
        d_prod = diff(prod_series)*100;
        d_prem = diff(wage_premium)*100;
        d_hours = diff(hours_series)*100;
        
        d_price = d_price+0.01*mean(price_series)*randn(size(diff(price_series)));
        d_price = -d_price;
        d_prod = d_prod+0.01*mean(prod_series)*randn(size(diff(prod_series)));
        d_prem = d_prem+0.01*mean(wage_premium)*randn(size(diff(wage_premium)));
        d_hours = d_hours+0.01*mean(hours_series)*randn(size(diff(hours_series)));
       
        if do_trend == 1
            [~,~,prod_trend] = bpass(d_prod,52,10e15);
            d_prod = d_prod-prod_trend;
            [~,~,hours_trend] = bpass(d_hours,52,10e15);
            d_hours = d_hours-hours_trend;
            [~,~,prem_trend] = bpass(d_prem,52,10e15);
            d_prem = d_prem-prem_trend;
            [~,~,price_trend] = bpass(d_price,52,10e15);
            d_price = d_price-price_trend;
        end
       
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % specification

%         % overidentifying restrictions (double zero)
        do_ov = 0;
%         % saving the shock series for comparison
%         do_save_shocks = 0;
        break_choice = 0;

        % order of the variables, i.e. identification choice
        % spec = 5
        data = [d_price d_prod d_prem d_hours];
        num_vars_diff = 4;
        nshocks = 2;
        num_vars_lev = num_vars - num_vars_diff;

        % more parameters for estimation

        prior = 1;                  % type of prior, 1: Litterman, 2: Flat, 3:Kadiyala, Karlsson
        nlagsvar = 8;               % lags in VAR
        c = 1;                      % c=1 constant, c=0: no constant, c>=2: trend with c dummy pieces
        ndraws = 1000;              % number of draws
        T = size(data,1)-nlagsvar;  % last observation to use in the estimation
        phi = [0.2 .5 10^5];        % RATS manual default coefficients for tightness on own first lag
        % other lags and exogenous variables
        decay = 3;                  % harmonic decay of own lags, 1:linear decay

        breaks = [];
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % estimating the BVAR

        disp('Step I: Estimating the BVAR')
        [Bet,Sigma,B,S] = bvarest(data,nlagsvar,num_vars_diff,c,breaks,break_choice,ndraws,prior,T,phi,decay);
        % Bet are draws from posterior, B is the mean of the posterior
        % S is Sols, Sigma is Sols for each of the ndraws

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % structural identification

        disp('Struct. Identification')
        if do_ov == 1
            break
        else
            [Response,Atilde,Respres]=ident(Bet,Sigma,nlagsvar,ndraws,HORIZON_ident,c,...
                num_vars,num_vars_diff,nshocks);
        end

        % sorting the Response
        for j=1:(num_vars*size(Bet,1))
            Respsort(j,:,:) = sort(Response(j,:,:),3);
        end;

        % Bayesian one-standard error
        Respmed(:,:) = Respsort(:,1:HORIZON_ident,0.5*ndraws);
        Resplow(:,:) = Respsort(:,1:HORIZON_ident,0.16*ndraws);
        Resphigh(:,:) = Respsort(:,1:HORIZON_ident,0.84*ndraws);

        Respmedall(:,:,it) = Respmed;
        Resplowall(:,:,it) = Resplow;
        Resphighall(:,:,it) = Resphigh;
    end % end iteration over simulation

    if do_trend == 1
        if prodfnc == 1
            save response_067 Respmedall
        elseif prodfnc == 2
            save response_117 Respmedall
        elseif prodfnc == 3
            save response_167 Respmedall
        elseif prodfnc == 4
            save response_217 Respmedall
        elseif prodfnc == 5
            save response_267 Respmedall
        elseif prodfnc == 6
            save response_317 Respmedall
        elseif prodfnc == 7
            save response_500 Respmedall
        end
    else
        if prodfnc == 1
            save response_067_not Respmedall
        elseif prodfnc == 2
            save response_117_not Respmedall
        elseif prodfnc == 3
            save response_167_not Respmedall
        elseif prodfnc == 4
            save response_217_not Respmedall
        elseif prodfnc == 5
            save response_267_not Respmedall
        elseif prodfnc == 6
            save response_317_not Respmedall
        elseif prodfnc == 7
            save response_500_not Respmedall
        end
    end
end
toc

for h=1:size(Respmedall,2)
    Respmed(:,h) = sum(Respmedall(:,h,:),3)/simsize;
    for n=1:size(Respmedall,1)
        Respstd(n,h) = std(Respmedall(n,h,:))/sqrt(simsize);
        Resplow(n,h) = Respmed(n,h)-Respstd(n,h);
        Resphigh(n,h) = Respmed(n,h)+Respstd(n,h);
    end
end
